####################################
# Code to Replicate Charts Figures and Tables
# Trust in Public Policy Algorithms
# Journal of Politics
# Geopolitical Experiments
# Study 1
# Ryan Kennedy, Philip D. Waggoner, and Matthew Ward
# 2021
####################################

# Load Needed Libraries
library(plyr)
library(arm)
library(ggplot2)
library(tidyr)
library(tibble)
library(gridExtra)
library(tidyverse)
library(stargazer)
library(grid)
library(here)

# Load the data from the first set of experiments
# Load the data
experiment1 <- read.csv(here("GeopoliticsExperiment1.csv"), as.is = TRUE)
experiment2 <- read.csv(here("GeopoliticsExperiment2.csv"), as.is = TRUE)
experiment3 <- read.csv(here("GeopoliticsExperiment3.csv"), as.is = TRUE)
experiment4 <- read.csv(here("GeopoliticsExperiment4.csv"), as.is = TRUE)

# Add in the labels for the experiment
experiment1$scenario_num <- 1
experiment2$scenario_num <- 2
experiment3$scenario_num <- 3
experiment4$scenario_num <- 4

# Values of advice for Experiments
experiment1$advice1 <- 1
experiment1$advice2 <- 9
experiment1$advice3 <- 50
experiment1$advice4 <- 28
experiment1$advice5 <- 2

experiment2$advice1 <- 79.5
experiment2$advice2 <- 20.5

experiment3$advice1 <- 25.6
experiment3$advice2 <- 25.7
experiment3$advice3 <- 34.0
experiment3$advice4 <- 14.7

experiment4$advice1 <- 11.6
experiment4$advice2 <- 60.7
experiment4$advice3 <- 13.4
experiment4$advice4 <- 14.3

### FIRST - Calculate Distance to Advice: abs(advice - userguess)
# Exp 1
experiment1$avgHuman.DISTANCE <- (abs(experiment1$advice1 - experiment1$Average.Human_1) +
                                    abs(experiment1$advice2 - experiment1$Average.Human_2) +
                                    abs(experiment1$advice3 - experiment1$Average.Human_3) +
                                    abs(experiment1$advice4 - experiment1$Average.Human_4) +
                                    abs(experiment1$advice5 - experiment1$Average.Human_5))/5

experiment1$randHuman.DISTANCE <- (abs(experiment1$advice1 - experiment1$Random.Human_1) +
                                     abs(experiment1$advice2 - experiment1$Random.Human_2) +
                                     abs(experiment1$advice3 - experiment1$Random.Human_3) +
                                     abs(experiment1$advice4 - experiment1$Random.Human_4) +
                                     abs(experiment1$advice5 - experiment1$Random.Human_5))/5

experiment1$algHuman.DISTANCE <- (abs(experiment1$advice1 - experiment1$Averaging.Algorithm_1) +
                                    abs(experiment1$advice2 - experiment1$Averaging.Algorithm_2) +
                                    abs(experiment1$advice3 - experiment1$Averaging.Algorithm_3) +
                                    abs(experiment1$advice4 - experiment1$Averaging.Algorithm_4) +
                                    abs(experiment1$advice5 - experiment1$Averaging.Algorithm_5))/5

experiment1$lerner.DISTANCE <- (abs(experiment1$advice1 - experiment1$Learning.Algorithm_1) +
                                  abs(experiment1$advice2 - experiment1$Learning.Algorithm_2) +
                                  abs(experiment1$advice3 - experiment1$Learning.Algorithm_3) +
                                  abs(experiment1$advice4 - experiment1$Learning.Algorithm_4) +
                                  abs(experiment1$advice5 - experiment1$Learning.Algorithm_5))/5

experiment1$distance <- rowMeans(cbind(experiment1$avgHuman.DISTANCE,experiment1$randHuman.DISTANCE,
                                       experiment1$algHuman.DISTANCE,experiment1$lerner.DISTANCE),na.rm=TRUE)

# Exp 2
experiment2$avgHuman.DISTANCE <- (abs(experiment2$advice1 - experiment2$Average.Human_1) +
                                    abs(experiment2$advice2 - experiment2$Average.Human_2))/2

experiment2$randHuman.DISTANCE <- (abs(experiment2$advice1 - experiment2$Random.Human_1) +
                                     abs(experiment2$advice2 - experiment2$Random.Human_2))/2

experiment2$algHuman.DISTANCE <- (abs(experiment2$advice1 - experiment2$Averaging.Algorithm_1) +
                                    abs(experiment2$advice2 - experiment2$Averaging.Algorithm_2))/2

experiment2$lerner.DISTANCE <- (abs(experiment2$advice1 - experiment2$Learning.Algorithm_1) +
                                  abs(experiment2$advice2 - experiment2$Learning.Algorithm_2))/2

experiment2$avgHuman.DISTANCE2 <- (abs(experiment2$advice1 - experiment2$Anch.Avg.Human_1) +
                                     abs(experiment2$advice2 - experiment2$Anch.Avg.Human_2))/2

experiment2$randHuman.DISTANCE2 <- (abs(experiment2$advice1 - experiment2$Anch.Rand.Human_1) +
                                      abs(experiment2$advice2 - experiment2$Anch.Rand.Human_2))/2

experiment2$algHuman.DISTANCE2 <- (abs(experiment2$advice1 - experiment2$Anch.alg.Combine_1) +
                                     abs(experiment2$advice2 - experiment2$Anch.alg.Combine_2))/2

experiment2$lerner.DISTANCE2 <- (abs(experiment2$advice1 - experiment2$Anch.Full.Alg_1) +
                                   abs(experiment2$advice2 - experiment2$Anch.Full.Alg_2))/2

experiment2$distance <- rowMeans(cbind(experiment2$avgHuman.DISTANCE,experiment2$randHuman.DISTANCE,
                                       experiment2$algHuman.DISTANCE,experiment2$lerner.DISTANCE,
                                       experiment2$avgHuman.DISTANCE2,experiment2$randHuman.DISTANCE2,
                                       experiment2$algHuman.DISTANCE2,experiment2$lerner.DISTANCE2),na.rm=TRUE)

# Exp 3
experiment3$avgHuman.DISTANCE <- (abs(experiment3$advice1 - experiment3$Average.Human_1) +
                                    abs(experiment3$advice2 - experiment3$Average.Human_2) +
                                    abs(experiment3$advice3 - experiment3$Average.Human_3) +
                                    abs(experiment3$advice4 - experiment3$Average.Human_4))/4

experiment3$randHuman.DISTANCE <- (abs(experiment3$advice1 - experiment3$Random.Human_1) +
                                     abs(experiment3$advice2 - experiment3$Random.Human_2) +
                                     abs(experiment3$advice3 - experiment3$Random.Human_3) +
                                     abs(experiment3$advice4 - experiment3$Random.Human_4))/4

experiment3$algHuman.DISTANCE <- (abs(experiment3$advice1 - experiment3$Averaging.Algorithm_1) +
                                    abs(experiment3$advice2 - experiment3$Averaging.Algorithm_2) +
                                    abs(experiment3$advice3 - experiment3$Averaging.Algorithm_3) +
                                    abs(experiment3$advice4 - experiment3$Averaging.Algorithm_4))/4

experiment3$lerner.DISTANCE <- (abs(experiment3$advice1 - experiment3$Learning.Algorithm_1) +
                                  abs(experiment3$advice2 - experiment3$Learning.Algorithm_2) +
                                  abs(experiment3$advice3 - experiment3$Learning.Algorithm_3) +
                                  abs(experiment3$advice4 - experiment3$Learning.Algorithm_4))/4

experiment3$avgHuman.DISTANCE2 <- (abs(experiment3$advice1 - experiment3$Anch.Avg.Human_1) +
                                     abs(experiment3$advice2 - experiment3$Anch.Avg.Human_2) +
                                     abs(experiment3$advice3 - experiment3$Anch.Avg.Human_3) +
                                     abs(experiment3$advice4 - experiment3$Anch.Avg.Human_4))/4

experiment3$randHuman.DISTANCE2 <- (abs(experiment3$advice1 - experiment3$Anch.Rand.Human_1) +
                                      abs(experiment3$advice2 - experiment3$Anch.Rand.Human_2) +
                                      abs(experiment3$advice3 - experiment3$Anch.Rand.Human_3) +
                                      abs(experiment3$advice4 - experiment3$Anch.Rand.Human_4))/4

experiment3$algHuman.DISTANCE2 <- (abs(experiment3$advice1 - experiment3$Anch.alg.Combine_1) +
                                     abs(experiment3$advice2 - experiment3$Anch.alg.Combine_2) +
                                     abs(experiment3$advice3 - experiment3$Anch.alg.Combine_3) +
                                     abs(experiment3$advice4 - experiment3$Anch.alg.Combine_4))/4

experiment3$lerner.DISTANCE2 <- (abs(experiment3$advice1 - experiment3$Anch.Full.Alg_1) +
                                   abs(experiment3$advice2 - experiment3$Anch.Full.Alg_2) +
                                   abs(experiment3$advice3 - experiment3$Anch.Full.Alg_3) +
                                   abs(experiment3$advice4 - experiment3$Anch.Full.Alg_4))/4

experiment3$distance <- rowMeans(cbind(experiment3$avgHuman.DISTANCE,experiment3$randHuman.DISTANCE,
                                       experiment3$algHuman.DISTANCE,experiment3$lerner.DISTANCE,
                                       experiment3$avgHuman.DISTANCE2,experiment3$randHuman.DISTANCE2,
                                       experiment3$algHuman.DISTANCE2,experiment3$lerner.DISTANCE2),na.rm=TRUE)


# Exp 4
experiment4$avgHuman.DISTANCE <- (abs(experiment4$advice1 - experiment4$Average.Human_1) +
                                    abs(experiment4$advice2 - experiment4$Average.Human_2) +
                                    abs(experiment4$advice3 - experiment4$Average.Human_3) +
                                    abs(experiment4$advice4 - experiment4$Average.Human_4))/4

experiment4$randHuman.DISTANCE <- (abs(experiment4$advice1 - experiment4$Random.Human_1) +
                                     abs(experiment4$advice2 - experiment4$Random.Human_2) +
                                     abs(experiment4$advice3 - experiment4$Random.Human_3) +
                                     abs(experiment4$advice4 - experiment4$Random.Human_4))/4

experiment4$algHuman.DISTANCE <- (abs(experiment4$advice1 - experiment4$Averaging.Algorithm_1) +
                                    abs(experiment4$advice2 - experiment4$Averaging.Algorithm_2) +
                                    abs(experiment4$advice3 - experiment4$Averaging.Algorithm_3) +
                                    abs(experiment4$advice4 - experiment4$Averaging.Algorithm_4))/4

experiment4$lerner.DISTANCE <- (abs(experiment4$advice1 - experiment4$Learning.Algorithm_1) +
                                  abs(experiment4$advice2 - experiment4$Learning.Algorithm_2) +
                                  abs(experiment4$advice3 - experiment4$Learning.Algorithm_3) +
                                  abs(experiment4$advice4 - experiment4$Learning.Algorithm_4))/4

experiment4$avgHuman.DISTANCE2 <- (abs(experiment4$advice1 - experiment4$Anch.Avg.Human_1) +
                                     abs(experiment4$advice2 - experiment4$Anch.Avg.Human_2) +
                                     abs(experiment4$advice3 - experiment4$Anch.Avg.Human_3) +
                                     abs(experiment4$advice4 - experiment4$Anch.Avg.Human_4))/4

experiment4$randHuman.DISTANCE2 <- (abs(experiment4$advice1 - experiment4$Anch.Rand.Human_1) +
                                      abs(experiment4$advice2 - experiment4$Anch.Rand.Human_2) +
                                      abs(experiment4$advice3 - experiment4$Anch.Rand.Human_3) +
                                      abs(experiment4$advice4 - experiment4$Anch.Rand.Human_4))/4

experiment4$algHuman.DISTANCE2 <- (abs(experiment4$advice1 - experiment4$Anch.alg.Combine_1) +
                                     abs(experiment4$advice2 - experiment4$Anch.alg.Combine_2) +
                                     abs(experiment4$advice3 - experiment4$Anch.alg.Combine_3) +
                                     abs(experiment4$advice4 - experiment4$Anch.alg.Combine_4))/4

experiment4$lerner.DISTANCE2 <- (abs(experiment4$advice1 - experiment4$Anch.Full.Alg_1) +
                                   abs(experiment4$advice2 - experiment4$Anch.Full.Alg_2) +
                                   abs(experiment4$advice3 - experiment4$Anch.Full.Alg_3) +
                                   abs(experiment4$advice4 - experiment4$Anch.Full.Alg_4))/4

experiment4$distance <- rowMeans(cbind(experiment4$avgHuman.DISTANCE,experiment4$randHuman.DISTANCE,
                                       experiment4$algHuman.DISTANCE,experiment4$lerner.DISTANCE,
                                       experiment4$avgHuman.DISTANCE2,experiment4$randHuman.DISTANCE2,
                                       experiment4$algHuman.DISTANCE2,experiment4$lerner.DISTANCE2),na.rm=TRUE)



### For second DV Calculate Weight of Advice: abs(userguesst2 - userguesst1)/abs(advice - userguesst1)
# Exp 1
experiment1$avgHumanWeight <- (abs(experiment1$Initial.Forecast_1 - experiment1$Average.Human_1) / abs(experiment1$advice1 - experiment1$Initial.Forecast_1) +
                                 abs(experiment1$Initial.Forecast_2 - experiment1$Average.Human_2) / abs(experiment1$advice2 - experiment1$Initial.Forecast_2) +
                                 abs(experiment1$Initial.Forecast_3 - experiment1$Average.Human_3) / abs(experiment1$advice3 - experiment1$Initial.Forecast_3) +
                                 abs(experiment1$Initial.Forecast_4 - experiment1$Average.Human_4) / abs(experiment1$advice4 - experiment1$Initial.Forecast_4) +
                                 abs(experiment1$Initial.Forecast_5 - experiment1$Average.Human_5) / abs(experiment1$advice5 - experiment1$Initial.Forecast_5))/5

experiment1$randHumanWeight <- (abs(experiment1$Initial.Forecast_1 - experiment1$Random.Human_1) / abs(experiment1$advice1 - experiment1$Initial.Forecast_1) +
                                  abs(experiment1$Initial.Forecast_2 - experiment1$Random.Human_2) / abs(experiment1$advice2 - experiment1$Initial.Forecast_2) +
                                  abs(experiment1$Initial.Forecast_3 - experiment1$Random.Human_3) / abs(experiment1$advice3 - experiment1$Initial.Forecast_3) +
                                  abs(experiment1$Initial.Forecast_4 - experiment1$Random.Human_4) / abs(experiment1$advice4 - experiment1$Initial.Forecast_4) +
                                  abs(experiment1$Initial.Forecast_5 - experiment1$Random.Human_5) / abs(experiment1$advice5 - experiment1$Initial.Forecast_5))/5

experiment1$algHumanWeight <- (abs(experiment1$Initial.Forecast_1 - experiment1$Averaging.Algorithm_1) / abs(experiment1$advice1 - experiment1$Initial.Forecast_1) +
                                 abs(experiment1$Initial.Forecast_2 - experiment1$Averaging.Algorithm_2) / abs(experiment1$advice2 - experiment1$Initial.Forecast_2) +
                                 abs(experiment1$Initial.Forecast_3 - experiment1$Averaging.Algorithm_3) / abs(experiment1$advice3 - experiment1$Initial.Forecast_3) +
                                 abs(experiment1$Initial.Forecast_4 - experiment1$Averaging.Algorithm_4) / abs(experiment1$advice4 - experiment1$Initial.Forecast_4) +
                                 abs(experiment1$Initial.Forecast_5 - experiment1$Averaging.Algorithm_5) / abs(experiment1$advice5 - experiment1$Initial.Forecast_5))/5

experiment1$lernerWeight <- (abs(experiment1$Initial.Forecast_1 - experiment1$Learning.Algorithm_1) / abs(experiment1$advice1 - experiment1$Initial.Forecast_1) +
                               abs(experiment1$Initial.Forecast_2 - experiment1$Learning.Algorithm_2) / abs(experiment1$advice2 - experiment1$Initial.Forecast_2) +
                               abs(experiment1$Initial.Forecast_3 - experiment1$Learning.Algorithm_3) / abs(experiment1$advice3 - experiment1$Initial.Forecast_3) +
                               abs(experiment1$Initial.Forecast_4 - experiment1$Learning.Algorithm_4) / abs(experiment1$advice4 - experiment1$Initial.Forecast_4) +
                               abs(experiment1$Initial.Forecast_5 - experiment1$Learning.Algorithm_5) / abs(experiment1$advice5 - experiment1$Initial.Forecast_5))/5

experiment1$adviceWt <- rowMeans(cbind(experiment1$avgHumanWeight,experiment1$randHumanWeight,
                                       experiment1$algHumanWeight,experiment1$lernerWeight),na.rm=TRUE)

# Exp 2
experiment2$avgHumanWeight <- (abs(experiment2$Initial.Forecast_1 - experiment2$Average.Human_1) / abs(experiment2$advice1 - experiment2$Initial.Forecast_1) +
                                 abs(experiment2$Initial.Forecast_2 - experiment2$Average.Human_2) / abs(experiment2$advice2 - experiment2$Initial.Forecast_2))/2

experiment2$randHumanWeight <- (abs(experiment2$Initial.Forecast_1 - experiment2$Random.Human_1) / abs(experiment2$advice1 - experiment2$Initial.Forecast_1) +
                                  abs(experiment2$Initial.Forecast_2 - experiment2$Random.Human_2) / abs(experiment2$advice2 - experiment2$Initial.Forecast_2))/2

experiment2$algHumanWeight <- (abs(experiment2$Initial.Forecast_1 - experiment2$Averaging.Algorithm_1) / abs(experiment2$advice1 - experiment2$Initial.Forecast_1) +
                                 abs(experiment2$Initial.Forecast_2 - experiment2$Averaging.Algorithm_2) / abs(experiment2$advice2 - experiment2$Initial.Forecast_2))/2

experiment2$lernerWeight <- (abs(experiment2$Initial.Forecast_1 - experiment2$Learning.Algorithm_1) / abs(experiment2$advice1 - experiment2$Initial.Forecast_1) +
                               abs(experiment2$Initial.Forecast_2 - experiment2$Learning.Algorithm_1) / abs(experiment2$advice2 - experiment2$Initial.Forecast_2))/2

experiment2$adviceWt <- rowMeans(cbind(experiment2$avgHumanWeight,experiment2$randHumanWeight,
                                       experiment2$algHumanWeight,experiment2$lernerWeight),na.rm=TRUE)


# Exp 3
experiment3$avgHumanWeight <- (abs(experiment3$Initial.Forecast_1 - experiment3$Average.Human_1) / abs(experiment3$advice1 - experiment3$Initial.Forecast_1) +
                                 abs(experiment3$Initial.Forecast_2 - experiment3$Average.Human_2) / abs(experiment3$advice2 - experiment3$Initial.Forecast_2) +
                                 abs(experiment3$Initial.Forecast_3 - experiment3$Average.Human_3) / abs(experiment3$advice3 - experiment3$Initial.Forecast_3) +
                                 abs(experiment3$Initial.Forecast_4 - experiment3$Average.Human_4) / abs(experiment3$advice4 - experiment3$Initial.Forecast_4))/4

experiment3$randHumanWeight <- (abs(experiment3$Initial.Forecast_1 - experiment3$Random.Human_1) / abs(experiment3$advice1 - experiment3$Initial.Forecast_1) +
                                  abs(experiment3$Initial.Forecast_2 - experiment3$Random.Human_2) / abs(experiment3$advice2 - experiment3$Initial.Forecast_2) +
                                  abs(experiment3$Initial.Forecast_3 - experiment3$Random.Human_3) / abs(experiment3$advice3 - experiment3$Initial.Forecast_3) +
                                  abs(experiment3$Initial.Forecast_4 - experiment3$Random.Human_4) / abs(experiment3$advice4 - experiment3$Initial.Forecast_4))/4

experiment3$algHumanWeight <- (abs(experiment3$Initial.Forecast_1 - experiment3$Averaging.Algorithm_1) / abs(experiment3$advice1 - experiment3$Initial.Forecast_1) +
                                 abs(experiment3$Initial.Forecast_2 - experiment3$Averaging.Algorithm_2) / abs(experiment3$advice2 - experiment3$Initial.Forecast_2) +
                                 abs(experiment3$Initial.Forecast_3 - experiment3$Averaging.Algorithm_3) / abs(experiment3$advice3 - experiment3$Initial.Forecast_3) +
                                 abs(experiment3$Initial.Forecast_4 - experiment3$Averaging.Algorithm_4) / abs(experiment3$advice4 - experiment3$Initial.Forecast_4))/4

experiment3$lernerWeight <- (abs(experiment3$Initial.Forecast_1 - experiment3$Learning.Algorithm_1) / abs(experiment3$advice1 - experiment3$Initial.Forecast_1) +
                               abs(experiment3$Initial.Forecast_2 - experiment3$Learning.Algorithm_2) / abs(experiment3$advice2 - experiment3$Initial.Forecast_2) +
                               abs(experiment3$Initial.Forecast_3 - experiment3$Learning.Algorithm_3) / abs(experiment3$advice3 - experiment3$Initial.Forecast_3) +
                               abs(experiment3$Initial.Forecast_4 - experiment3$Learning.Algorithm_4) / abs(experiment3$advice4 - experiment3$Initial.Forecast_4))/4

experiment3$adviceWt <- rowMeans(cbind(experiment3$avgHumanWeight,experiment3$randHumanWeight,
                                       experiment3$algHumanWeight,experiment3$lernerWeight),na.rm=TRUE)

# Exp 4
experiment4$avgHumanWeight <- (abs(experiment4$Initial.Forecast_1 - experiment4$Average.Human_1) / abs(experiment4$advice1 - experiment4$Initial.Forecast_1) +
                                 abs(experiment4$Initial.Forecast_2 - experiment4$Average.Human_2) / abs(experiment4$advice2 - experiment4$Initial.Forecast_2) +
                                 abs(experiment4$Initial.Forecast_3 - experiment4$Average.Human_3) / abs(experiment4$advice3 - experiment4$Initial.Forecast_3) +
                                 abs(experiment4$Initial.Forecast_4 - experiment4$Average.Human_4) / abs(experiment4$advice4 - experiment4$Initial.Forecast_4))/4

experiment4$randHumanWeight <- (abs(experiment4$Initial.Forecast_1 - experiment4$Random.Human_1) / abs(experiment4$advice1 - experiment4$Initial.Forecast_1) +
                                  abs(experiment4$Initial.Forecast_2 - experiment4$Random.Human_2) / abs(experiment4$advice2 - experiment4$Initial.Forecast_2) +
                                  abs(experiment4$Initial.Forecast_3 - experiment4$Random.Human_3) / abs(experiment4$advice3 - experiment4$Initial.Forecast_3) +
                                  abs(experiment4$Initial.Forecast_4 - experiment4$Random.Human_4) / abs(experiment4$advice4 - experiment4$Initial.Forecast_4))/4

experiment4$algHumanWeight <- (abs(experiment4$Initial.Forecast_1 - experiment4$Averaging.Algorithm_1) / abs(experiment4$advice1 - experiment4$Initial.Forecast_1) +
                                 abs(experiment4$Initial.Forecast_2 - experiment4$Averaging.Algorithm_2) / abs(experiment4$advice2 - experiment4$Initial.Forecast_2) +
                                 abs(experiment4$Initial.Forecast_3 - experiment4$Averaging.Algorithm_3) / abs(experiment4$advice3 - experiment4$Initial.Forecast_3) +
                                 abs(experiment4$Initial.Forecast_4 - experiment4$Averaging.Algorithm_4) / abs(experiment4$advice4 - experiment4$Initial.Forecast_4))/4

experiment4$lernerWeight <- (abs(experiment4$Initial.Forecast_1 - experiment4$Learning.Algorithm_1) / abs(experiment4$advice1 - experiment4$Initial.Forecast_1) +
                               abs(experiment4$Initial.Forecast_2 - experiment4$Learning.Algorithm_2) / abs(experiment4$advice2 - experiment4$Initial.Forecast_2) +
                               abs(experiment4$Initial.Forecast_3 - experiment4$Learning.Algorithm_3) / abs(experiment4$advice3 - experiment4$Initial.Forecast_3) +
                               abs(experiment4$Initial.Forecast_4 - experiment4$Learning.Algorithm_4) / abs(experiment4$advice4 - experiment4$Initial.Forecast_4))/4

experiment4$adviceWt <- rowMeans(cbind(experiment4$avgHumanWeight,experiment4$randHumanWeight,
                                       experiment4$algHumanWeight,experiment4$lernerWeight),na.rm=TRUE)


##  Calculate Brier scores as robustness check/secondary analysis
# Calculate Brier Distance to Advice
# Exp 1
experiment1$avgHumanBrier <- ((0.01*experiment1$Average.Human_1 - 0.01*experiment1$advice1)^2 +
                                (0.01*experiment1$Average.Human_2 - 0.01*experiment1$advice2)^2 +
                                (0.01*experiment1$Average.Human_3 - 0.01*experiment1$advice3)^2 +
                                (0.01*experiment1$Average.Human_4 - 0.01*experiment1$advice4)^2 +
                                (0.01*experiment1$Average.Human_5 - 0.01*experiment1$advice5)^2)/5

experiment1$randHumanBrier <- ((0.01*experiment1$Random.Human_1 - 0.01*experiment1$advice1)^2 +
                                 (0.01*experiment1$Random.Human_2 - 0.01*experiment1$advice2)^2 +
                                 (0.01*experiment1$Random.Human_3 - 0.01*experiment1$advice3)^2 +
                                 (0.01*experiment1$Random.Human_4 - 0.01*experiment1$advice4)^2 +
                                 (0.01*experiment1$Random.Human_5 - 0.01*experiment1$advice5)^2)/5

experiment1$algHumanBrier <- ((0.01*experiment1$Averaging.Algorithm_1 - 0.01*experiment1$advice1)^2 +
                                (0.01*experiment1$Averaging.Algorithm_2 - 0.01*experiment1$advice2)^2 +
                                (0.01*experiment1$Averaging.Algorithm_3 - 0.01*experiment1$advice3)^2 +
                                (0.01*experiment1$Averaging.Algorithm_4 - 0.01*experiment1$advice4)^2 +
                                (0.01*experiment1$Averaging.Algorithm_5 - 0.01*experiment1$advice5)^2)/5

experiment1$lernerBrier <- ((0.01*experiment1$Learning.Algorithm_1 - 0.01*experiment1$advice1)^2 +
                              (0.01*experiment1$Learning.Algorithm_2 - 0.01*experiment1$advice2)^2 +
                              (0.01*experiment1$Learning.Algorithm_3 - 0.01*experiment1$advice3)^2 +
                              (0.01*experiment1$Learning.Algorithm_4 - 0.01*experiment1$advice4)^2 +
                              (0.01*experiment1$Learning.Algorithm_5 - 0.01*experiment1$advice5)^2)/5

experiment1$brier <- rowMeans(cbind(experiment1$avgHumanBrier,experiment1$randHumanBrier,
                                    experiment1$algHumanBrier,experiment1$lernerBrier),na.rm=TRUE)

# Exp 2
experiment2$avgHumanBrier <- ((0.01*experiment2$Average.Human_1 - 0.01*experiment2$advice1)^2 +
                                (0.01*experiment2$Average.Human_2 - 0.01*experiment2$advice2)^2)/2

experiment2$randHumanBrier <- ((0.01*experiment2$Random.Human_1 - 0.01*experiment2$advice1)^2 +
                                 (0.01*experiment2$Random.Human_2 - 0.01*experiment2$advice2)^2)/2

experiment2$algHumanBrier <- ((0.01*experiment2$Averaging.Algorithm_1 - 0.01*experiment2$advice1)^2 +
                                (0.01*experiment2$Averaging.Algorithm_2 - 0.01*experiment2$advice2)^2)/2

experiment2$lernerBrier <- ((0.01*experiment2$Learning.Algorithm_1 - 0.01*experiment2$advice1)^2 +
                              (0.01*experiment2$Learning.Algorithm_2 - 0.01*experiment2$advice2)^2)/2

experiment2$avgHumanBrier2 <- ((0.01*experiment2$Anch.Avg.Human_1 - 0.01*experiment2$advice1)^2 +
                                 (0.01*experiment2$Anch.Avg.Human_2 - 0.01*experiment2$advice2)^2)/2

experiment2$randHumanBrier2 <- ((0.01*experiment2$Anch.Rand.Human_1 - 0.01*experiment2$advice1)^2 +
                                  (0.01*experiment2$Anch.Rand.Human_2 - 0.01*experiment2$advice2)^2)/2

experiment2$algHumanBrier2 <- ((0.01*experiment2$Anch.alg.Combine_1 - 0.01*experiment2$advice1)^2 +
                                 (0.01*experiment2$Anch.alg.Combine_2 - 0.01*experiment2$advice2)^2)/2

experiment2$lernerBrier2 <- ((0.01*experiment2$Anch.Full.Alg_1 - 0.01*experiment2$advice1)^2 +
                               (0.01*experiment2$Anch.Full.Alg_2 - 0.01*experiment2$advice2)^2)/2

experiment2$brier <- rowMeans(cbind(experiment2$avgHumanBrier,experiment2$randHumanBrier,
                                    experiment2$algHumanBrier,experiment2$lernerBrier,
                                    experiment2$avgHumanBrier2,experiment2$randHumanBrier2,
                                    experiment2$algHumanBrier2,experiment2$lernerBrier2),na.rm=TRUE)

# Exp 3
experiment3$avgHumanBrier <- ((0.01*experiment3$Average.Human_1 - 0.01*experiment3$advice1)^2 +
                                (0.01*experiment3$Average.Human_2 - 0.01*experiment3$advice2)^2 +
                                (0.01*experiment3$Average.Human_3 - 0.01*experiment3$advice3)^2 +
                                (0.01*experiment3$Average.Human_4 - 0.01*experiment3$advice4)^2)/4

experiment3$randHumanBrier <- ((0.01*experiment3$Random.Human_1 - 0.01*experiment3$advice1)^2 +
                                 (0.01*experiment3$Random.Human_2 - 0.01*experiment3$advice2)^2 +
                                 (0.01*experiment3$Random.Human_3 - 0.01*experiment3$advice3)^2 +
                                 (0.01*experiment3$Random.Human_4 - 0.01*experiment3$advice4)^2)/4

experiment3$algHumanBrier <- ((0.01*experiment3$Averaging.Algorithm_1 - 0.01*experiment3$advice1)^2 +
                                (0.01*experiment3$Averaging.Algorithm_2 - 0.01*experiment3$advice2)^2 +
                                (0.01*experiment3$Averaging.Algorithm_3 - 0.01*experiment3$advice3)^2 +
                                (0.01*experiment3$Averaging.Algorithm_4 - 0.01*experiment3$advice4)^2)/4

experiment3$lernerBrier <- ((0.01*experiment3$Learning.Algorithm_1 - 0.01*experiment3$advice1)^2 +
                              (0.01*experiment3$Learning.Algorithm_2 - 0.01*experiment3$advice2)^2 +
                              (0.01*experiment3$Learning.Algorithm_3 - 0.01*experiment3$advice3)^2 +
                              (0.01*experiment3$Learning.Algorithm_4 - 0.01*experiment3$advice4)^2)/4

experiment3$avgHumanBrier2 <- ((0.01*experiment3$Anch.Avg.Human_1 - 0.01*experiment3$advice1)^2 +
                                 (0.01*experiment3$Anch.Avg.Human_2 - 0.01*experiment3$advice2)^2 +
                                 (0.01*experiment3$Anch.Avg.Human_3 - 0.01*experiment3$advice3)^2 +
                                 (0.01*experiment3$Anch.Avg.Human_4 - 0.01*experiment3$advice4)^2)/4

experiment3$randHumanBrier2 <- ((0.01*experiment3$Anch.Rand.Human_1 - 0.01*experiment3$advice1)^2 +
                                  (0.01*experiment3$Anch.Rand.Human_2 - 0.01*experiment3$advice2)^2 +
                                  (0.01*experiment3$Anch.Rand.Human_3 - 0.01*experiment3$advice3)^2 +
                                  (0.01*experiment3$Anch.Rand.Human_4 - 0.01*experiment3$advice4)^2)/4

experiment3$algHumanBrier2 <- ((0.01*experiment3$Anch.alg.Combine_1 - 0.01*experiment3$advice1)^2 +
                                 (0.01*experiment3$Anch.alg.Combine_2 - 0.01*experiment3$advice2)^2 +
                                 (0.01*experiment3$Anch.alg.Combine_3 - 0.01*experiment3$advice3)^2 +
                                 (0.01*experiment3$Anch.alg.Combine_4 - 0.01*experiment3$advice4)^2)/4

experiment3$lernerBrier2 <- ((0.01*experiment3$Anch.Full.Alg_1 - 0.01*experiment3$advice1)^2 +
                               (0.01*experiment3$Anch.Full.Alg_2 - 0.01*experiment3$advice2)^2 +
                               (0.01*experiment3$Anch.Full.Alg_3 - 0.01*experiment3$advice3)^2 +
                               (0.01*experiment3$Anch.Full.Alg_4 - 0.01*experiment3$advice4)^2)/4

experiment3$brier <- rowMeans(cbind(experiment3$avgHumanBrier,experiment3$randHumanBrier,
                                    experiment3$algHumanBrier,experiment3$lernerBrier,
                                    experiment3$avgHumanBrier2,experiment3$randHumanBrier2,
                                    experiment3$algHumanBrier2,experiment3$lernerBrier2),na.rm=TRUE)

# Exp 4
experiment4$avgHumanBrier <- ((0.01*experiment4$Average.Human_1 - 0.01*experiment4$advice1)^2 +
                                (0.01*experiment4$Average.Human_2 - 0.01*experiment4$advice2)^2 +
                                (0.01*experiment4$Average.Human_3 - 0.01*experiment4$advice3)^2 +
                                (0.01*experiment4$Average.Human_4 - 0.01*experiment4$advice4)^2)/4

experiment4$randHumanBrier <- ((0.01*experiment4$Random.Human_1 - 0.01*experiment4$advice1)^2 +
                                 (0.01*experiment4$Random.Human_2 - 0.01*experiment4$advice2)^2 +
                                 (0.01*experiment4$Random.Human_3 - 0.01*experiment4$advice3)^2 +
                                 (0.01*experiment4$Random.Human_4 - 0.01*experiment4$advice4)^2)/4

experiment4$algHumanBrier <- ((0.01*experiment4$Averaging.Algorithm_1 - 0.01*experiment4$advice1)^2 +
                                (0.01*experiment4$Averaging.Algorithm_2 - 0.01*experiment4$advice2)^2 +
                                (0.01*experiment4$Averaging.Algorithm_3 - 0.01*experiment4$advice3)^2 +
                                (0.01*experiment4$Averaging.Algorithm_4 - 0.01*experiment4$advice4)^2)/4

experiment4$lernerBrier <- ((0.01*experiment4$Learning.Algorithm_1 - 0.01*experiment4$advice1)^2 +
                              (0.01*experiment4$Learning.Algorithm_2 - 0.01*experiment4$advice2)^2 +
                              (0.01*experiment4$Learning.Algorithm_3 - 0.01*experiment4$advice3)^2 +
                              (0.01*experiment4$Learning.Algorithm_4 - 0.01*experiment4$advice4)^2)/4

experiment4$avgHumanBrier2 <- ((0.01*experiment4$Anch.Avg.Human_1 - 0.01*experiment4$advice1)^2 +
                                 (0.01*experiment4$Anch.Avg.Human_2 - 0.01*experiment4$advice2)^2 +
                                 (0.01*experiment4$Anch.Avg.Human_3 - 0.01*experiment4$advice3)^2 +
                                 (0.01*experiment4$Anch.Avg.Human_4 - 0.01*experiment4$advice4)^2)/4

experiment4$randHumanBrier2 <- ((0.01*experiment4$Anch.Rand.Human_1 - 0.01*experiment4$advice1)^2 +
                                  (0.01*experiment4$Anch.Rand.Human_2 - 0.01*experiment4$advice2)^2 +
                                  (0.01*experiment4$Anch.Rand.Human_3 - 0.01*experiment4$advice3)^2 +
                                  (0.01*experiment4$Anch.Rand.Human_4 - 0.01*experiment4$advice4)^2)/4

experiment4$algHumanBrier2 <- ((0.01*experiment4$Anch.alg.Combine_1 - 0.01*experiment4$advice1)^2 +
                                 (0.01*experiment4$Anch.alg.Combine_2 - 0.01*experiment4$advice2)^2 +
                                 (0.01*experiment4$Anch.alg.Combine_3 - 0.01*experiment4$advice3)^2 +
                                 (0.01*experiment4$Anch.alg.Combine_4 - 0.01*experiment4$advice4)^2)/4

experiment4$lernerBrier2 <- ((0.01*experiment4$Anch.Full.Alg_1 - 0.01*experiment4$advice1)^2 +
                               (0.01*experiment4$Anch.Full.Alg_2 - 0.01*experiment4$advice2)^2 +
                               (0.01*experiment4$Anch.Full.Alg_3 - 0.01*experiment4$advice3)^2 +
                               (0.01*experiment4$Anch.Full.Alg_4 - 0.01*experiment4$advice4)^2)/4

experiment4$brier <- rowMeans(cbind(experiment4$avgHumanBrier,experiment4$randHumanBrier,
                                    experiment4$algHumanBrier,experiment4$lernerBrier,
                                    experiment4$avgHumanBrier2,experiment4$randHumanBrier2,
                                    experiment4$algHumanBrier2,experiment4$lernerBrier2),na.rm=TRUE)


# Combine the data from the experiments
experiments <- rbind.fill(experiment1, experiment2, experiment3, experiment4)


# Add in treatment indicators
experiments$AvgHumanTreat <- ifelse(!is.na(experiments$avgHuman.DISTANCE) | 
                                      !is.na(experiments$avgHuman.DISTANCE2) |
                                      !is.na(experiments$avgHumanWeight), 1, 0) 

experiments$RandHumanTreat <- ifelse(!is.na(experiments$randHuman.DISTANCE) | 
                                       !is.na(experiments$randHuman.DISTANCE2) |
                                       !is.na(experiments$randHumanWeight), 1, 0) 

experiments$AlgHumanTreat <- ifelse(!is.na(experiments$algHuman.DISTANCE) | 
                                      !is.na(experiments$algHuman.DISTANCE2) |
                                      !is.na(experiments$algHumanWeight), 1, 0) 

experiments$LernerTreat <- ifelse(!is.na(experiments$lerner.DISTANCE) | 
                                    !is.na(experiments$lerner.DISTANCE2) |
                                    !is.na(experiments$lernerWeight), 1, 0) 

experiments$Anchoring <- ifelse(!is.na(experiments$avgHuman.DISTANCE2) |
                                  !is.na(experiments$randHuman.DISTANCE2) |
                                  !is.na(experiments$algHuman.DISTANCE2) |
                                  !is.na(experiments$lerner.DISTANCE2) |
                                  !is.na(experiments$avgHumanBrier2) |
                                  !is.na(experiments$randHumanBrier2) |
                                  !is.na(experiments$algHumanBrier2) |
                                  !is.na(experiments$lernerBrier2), 1, 0)


# Gender
experiments$female[experiments$Gender == "Male"] <- 0
experiments$female[experiments$Gender == "Female"] <- 1

# Education
experiments$edvalue[experiments$education == "High school graduate/GED" ] <- 1
experiments$edvalue[experiments$education == "Trade or vocational certification" ] <- 1
experiments$edvalue[experiments$education == "Some college, or an associate degree"] <- 2
experiments$edvalue[experiments$education == "Bachelor's degree (for example, BA, AB, BS)"] <- 2
experiments$edvalue[experiments$education == "Some graduate school" ] <- 3
experiments$edvalue[experiments$education == "Master's degree (for example, MA, MSW, MBA)" ] <- 3
experiments$edvalue[experiments$education == "Professional degree (for example, MD, JD, DDS)" ] <- 3
experiments$edvalue[experiments$education == "Doctoral degree (for example, PhD, EdD)" ] <- 3
experiments$edvalue[experiments$scenario_num == 1 | experiments$scenario_num == 2] <- 2 #Student samples, so, by definition, value of 2

experiments$ed <- experiments$edvalue


# Trust in Automation
experiments$toa1[experiments$Trust.Automation_1 == "Strongly agree"] <- 7
experiments$toa1[experiments$Trust.Automation_1 == "Agree"] <- 6
experiments$toa1[experiments$Trust.Automation_1 == "Somewhat agree"] <- 5
experiments$toa1[experiments$Trust.Automation_1 == "Neither agree nor disagree"] <- 4
experiments$toa1[experiments$Trust.Automation_1 == "Somewhat disagree"] <- 3
experiments$toa1[experiments$Trust.Automation_1 == "Disagree"] <- 2
experiments$toa1[experiments$Trust.Automation_1 == "Strongly disagree"] <- 1

experiments$toa2[experiments$Trust.Automation_2 == "Strongly agree"] <- 7
experiments$toa2[experiments$Trust.Automation_2 == "Agree"] <- 6
experiments$toa2[experiments$Trust.Automation_2 == "Somewhat agree"] <- 5
experiments$toa2[experiments$Trust.Automation_2 == "Neither agree nor disagree"] <- 4
experiments$toa2[experiments$Trust.Automation_2 == "Somewhat disagree"] <- 3
experiments$toa2[experiments$Trust.Automation_2 == "Disagree"] <- 2
experiments$toa2[experiments$Trust.Automation_2 == "Strongly disagree"] <- 1

experiments$toa3[experiments$Trust.Automation_3 == "Strongly agree"] <- 7
experiments$toa3[experiments$Trust.Automation_3 == "Agree"] <- 6
experiments$toa3[experiments$Trust.Automation_3 == "Somewhat agree"] <- 5
experiments$toa3[experiments$Trust.Automation_3 == "Neither agree nor disagree"] <- 4
experiments$toa3[experiments$Trust.Automation_3 == "Somewhat disagree"] <- 3
experiments$toa3[experiments$Trust.Automation_3 == "Disagree"] <- 2
experiments$toa3[experiments$Trust.Automation_3 == "Strongly disagree"] <- 1

experiments$toa4[experiments$Trust.Automation_4 == "Strongly agree"] <- 1
experiments$toa4[experiments$Trust.Automation_4 == "Agree"] <- 2
experiments$toa4[experiments$Trust.Automation_4 == "Somewhat agree"] <- 3
experiments$toa4[experiments$Trust.Automation_4 == "Neither agree nor disagree"] <- 4
experiments$toa4[experiments$Trust.Automation_4 == "Somewhat disagree"] <- 5
experiments$toa4[experiments$Trust.Automation_4 == "Disagree"] <- 6
experiments$toa4[experiments$Trust.Automation_4 == "Strongly disagree"] <- 7

experiments$toa5[experiments$Trust.Automation_5 == "Strongly agree"] <- 1
experiments$toa5[experiments$Trust.Automation_5 == "Agree"] <- 2
experiments$toa5[experiments$Trust.Automation_5 == "Somewhat agree"] <- 3
experiments$toa5[experiments$Trust.Automation_5 == "Neither agree nor disagree"] <- 4
experiments$toa5[experiments$Trust.Automation_5 == "Somewhat disagree"] <- 5
experiments$toa5[experiments$Trust.Automation_5 == "Disagree"] <- 6
experiments$toa5[experiments$Trust.Automation_5 == "Strongly disagree"] <- 7

experiments$toa6[experiments$Trust.Automation_6 == "Strongly agree"] <- 7
experiments$toa6[experiments$Trust.Automation_6 == "Agree"] <- 6
experiments$toa6[experiments$Trust.Automation_6 == "Somewhat agree"] <- 5
experiments$toa6[experiments$Trust.Automation_6 == "Neither agree nor disagree"] <- 4
experiments$toa6[experiments$Trust.Automation_6 == "Somewhat disagree"] <- 3
experiments$toa6[experiments$Trust.Automation_6 == "Disagree"] <- 2
experiments$toa6[experiments$Trust.Automation_6 == "Strongly disagree"] <- 1

experiments$toa7[experiments$Trust.Automation_7 == "Strongly agree"] <- 7
experiments$toa7[experiments$Trust.Automation_7 == "Agree"] <- 6
experiments$toa7[experiments$Trust.Automation_7 == "Somewhat agree"] <- 5
experiments$toa7[experiments$Trust.Automation_7 == "Neither agree nor disagree"] <- 4
experiments$toa7[experiments$Trust.Automation_7 == "Somewhat disagree"] <- 3
experiments$toa7[experiments$Trust.Automation_7 == "Disagree"] <- 2
experiments$toa7[experiments$Trust.Automation_7 == "Strongly disagree"] <- 1

experiments$toa8[experiments$Trust.Automation_8 == "Strongly agree"] <- 1
experiments$toa8[experiments$Trust.Automation_8 == "Agree"] <- 2
experiments$toa8[experiments$Trust.Automation_8 == "Somewhat agree"] <- 3
experiments$toa8[experiments$Trust.Automation_8 == "Neither agree nor disagree"] <- 4
experiments$toa8[experiments$Trust.Automation_8 == "Somewhat disagree"] <- 5
experiments$toa8[experiments$Trust.Automation_8 == "Disagree"] <- 6
experiments$toa8[experiments$Trust.Automation_8 == "Strongly disagree"] <- 7

experiments$trustAutomation = (experiments$toa1 + experiments$toa2 + experiments$toa3 + experiments$toa4 + experiments$toa5 + experiments$toa6 + experiments$toa7 + experiments$toa8)/8

experiments$tia <- experiments$trustAutomation

# Age
experiments$age <- experiments$Age

#### Create a table of the summary statistics
# Table A1 in SI
experimentsSum <- experiments %>%
  select(distance, brier, RandHumanTreat, AvgHumanTreat,
         AlgHumanTreat, LernerTreat, Anchoring, tia, age, ed, female)

stargazer(experimentsSum)

# Note, the 2nd row in the table is inserted manually to reflect characteristics of weight of advice
experiments %>%
  filter(adviceWt <= 1) %>%
  select(adviceWt) %>%
  summarize(n = n(),
            mn = mean(adviceWt),
            stdev = sd(adviceWt),
            mnm = min(adviceWt),
            mx = max(adviceWt))


#### Create Balance Table 
# See Table A5 in SI
# Compare means for each experimental condition
experiments %>% 
  mutate(treatment = case_when(RandHumanTreat == 1 ~ "RH",
                               AvgHumanTreat == 1 ~ "AH",
                               AlgHumanTreat == 1 ~ "ApH",
                               LernerTreat == 1 ~ "A")) %>%
  filter(!is.na(treatment)) %>%
  group_by(treatment) %>%
  summarise(med = mean(ed, na.rm = T), 
            mage = mean(age, na.rm = T),
            mf = mean(female, na.rm = T)) %>%
  gather(key = "variable","value",med,mage,mf) %>%
  spread(treatment,value)

# T-tests for differences in means between experimental conditions
experiments %>%
  mutate(treatment = case_when(RandHumanTreat == 1 ~ "RH",
                               AvgHumanTreat == 1 ~ "AH",
                               AlgHumanTreat == 1 ~ "ApH",
                               LernerTreat == 1 ~ "A")) %>%
  filter(treatment == "RH" | treatment == "AH") %>%
  select(treatment, ed, age, female) %>%
  gather(key = variable, value = value, -treatment) %>%
  group_by(treatment, variable) %>%
  summarise(value = list(value)) %>%
  spread(treatment, value) %>%
  group_by(variable) %>%
  mutate(p_value = t.test(unlist(RH), unlist(AH))$p.value)

experiments %>%
  mutate(treatment = case_when(RandHumanTreat == 1 ~ "RH",
                               AvgHumanTreat == 1 ~ "AH",
                               AlgHumanTreat == 1 ~ "ApH",
                               LernerTreat == 1 ~ "A")) %>%
  filter(treatment == "RH" | treatment == "ApH") %>%
  select(treatment, ed, age, female) %>%
  gather(key = variable, value = value, -treatment) %>%
  group_by(treatment, variable) %>%
  summarise(value = list(value)) %>%
  spread(treatment, value) %>%
  group_by(variable) %>%
  mutate(p_value = t.test(unlist(RH), unlist(ApH))$p.value)

experiments %>%
  mutate(treatment = case_when(RandHumanTreat == 1 ~ "RH",
                               AvgHumanTreat == 1 ~ "AH",
                               AlgHumanTreat == 1 ~ "ApH",
                               LernerTreat == 1 ~ "A")) %>%
  filter(treatment == "RH" | treatment == "A") %>%
  select(treatment, ed, age, female) %>%
  gather(key = variable, value = value, -treatment) %>%
  group_by(treatment, variable) %>%
  summarise(value = list(value)) %>%
  spread(treatment, value) %>%
  group_by(variable) %>%
  mutate(p_value = t.test(unlist(RH), unlist(A))$p.value)

#### Descriptive analysis of advice weight
# Mean and median advice weight for Study 1
experiments %>%
  filter(adviceWt <= 1) %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 

experiments %>%
  filter(adviceWt <= 1 & scenario_num == 1) %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 

experiments %>%
  filter(adviceWt <= 1 & scenario_num == 2) %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 

experiments %>%
  filter(adviceWt <= 1 & scenario_num == 3) %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 

experiments %>%
  filter(adviceWt <= 1 & scenario_num == 4) %>%
  summarize(avg = mean(adviceWt, na.rm = T),
            med = median(adviceWt, na.rm = T)) 

# Figure A1 in SI
experiments %>%
  filter(adviceWt <= 1) %>%
  mutate(adType = case_when(RandHumanTreat == 1 ~ "Random Respondent",
                            AvgHumanTreat == 1 ~ "Average of Respondents",
                            AlgHumanTreat == 1 ~ "Algorithm Aggregation",
                            LernerTreat == 1 ~ "Algorithm")) %>%
  ggplot() + geom_histogram(aes(x = adviceWt, fill = adType), binwidth = .1) + 
  facet_wrap(~adType, ncol = 2) + theme_bw() +
  labs(x = "Advice Weight", y = "Count", fill = "Source")

# Figure A2 in SI
acc1 <- experiments %>%
  filter(adviceWt <= 1) %>%
  mutate(acceptAdvice = ifelse(adviceWt > 0.5, 1, 0),
         adType = case_when(RandHumanTreat == 1 ~ "Random Respondent",
                            AvgHumanTreat == 1 ~ "Average of Respondents",
                            AlgHumanTreat == 1 ~ "Algorithm Aggregation",
                            LernerTreat == 1 ~ "Algorithm"),
         Experiment = case_when(scenario_num == 1 ~ "Experiment 1",
                                scenario_num == 2 ~ "Experiment 2",
                                scenario_num == 3 ~ "Experiment 3",
                                scenario_num == 4 ~ "Experiment 4")) %>%
  group_by(adType, Experiment) %>%
  summarise(respondents = n(),
            sumAccept = sum(acceptAdvice)) %>%
  ungroup() %>%
  mutate(pctAccept = sumAccept/respondents) %>%
  ggplot() + geom_bar(aes(x = Experiment, y = pctAccept, fill = adType), stat = "identity", position = "dodge") +
  #scale_y_continuous(breaks = seq(0, 0.3, 0.1)) +
  labs(x = "Experiment", y = "Percent Accepting Advice (Weight > 0.5)", fill = "Source") +
  theme_bw()
acc1

#### MODELS
model1b <- lmer(distance ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + (1|scenario_num),
                data = experiments); summary(model1b)

model2b <- lmer(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + (1|scenario_num),
                data = experiments[experiments$adviceWt <= 1,]); summary(model2b)

model3b <- lmer(brier ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + (1|scenario_num),
                data = experiments); summary(model3b)

#### PLOTS
# Simulation for distance plot
# Figure 1
mod1sum <- summary(model1b)
mod1sim <- sim(model1b, n.sims = 1000)
mod1simcoef <- as_tibble(data.frame(coef(mod1sim)$fixef))
mod1simcoef$simulation <- rownames(mod1simcoef)
mod1simcoeflong <- gather(mod1simcoef, condition, coefficient, X.Intercept.:LernerTreat)
mod1simcoeflong$condition[mod1simcoeflong$condition == "LernerTreat"] <- "Computer Algorithm"
mod1simcoeflong$condition[mod1simcoeflong$condition == "AlgHumanTreat"] <- "Algorithm Aggregation"
mod1simcoeflong$condition[mod1simcoeflong$condition == "AvgHumanTreat"] <- "Average of Humans"

# Plot of Figure 1
m1sim95 <- mod1simcoeflong %>%
  filter(condition %in% c("Computer Algorithm", "Algorithm Aggregation", "Average of Humans")) %>%
  group_by(condition) %>%
  summarize(mcond = mean(coefficient, na.rm = T),
            sdcond = sd(coefficient, na.rm = T)) %>%
  mutate(hi = mcond + 1.96*sdcond,
         lo = mcond - 1.96*sdcond) %>% 
  ggplot() +
  geom_point(aes(x = condition, y = mcond), size = 3, color = "red") +
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(x = condition, ymin = lo, ymax = hi), color = "red") +
  labs(x = "Condition", y = "ATE", title = "Distance to Advice") +
  coord_flip() + theme_bw()
m1sim95


# Simulation for advice weight plot\
# Figure 1
mod2sum <- summary(model2b)
mod2sim <- sim(model2b, n.sims = 1000)
mod2simcoef <- as_tibble(data.frame(coef(mod2sim)$fixef))
mod2simcoef$simulation <- rownames(mod2simcoef)
mod2simcoeflong <- gather(mod2simcoef, condition, coefficient, X.Intercept.:LernerTreat)
mod2simcoeflong$condition[mod2simcoeflong$condition == "LernerTreat"] <- "Computer Algorithm"
mod2simcoeflong$condition[mod2simcoeflong$condition == "AlgHumanTreat"] <- "Algorithm Aggregation"
mod2simcoeflong$condition[mod2simcoeflong$condition == "AvgHumanTreat"] <- "Average of Humans"

m2sim95 <- mod2simcoeflong %>%
  filter(condition %in% c("Computer Algorithm", "Algorithm Aggregation", "Average of Humans")) %>%
  group_by(condition) %>%
  summarize(mcond = mean(coefficient, na.rm = T),
            sdcond = sd(coefficient, na.rm = T)) %>%
  mutate(hi = mcond + 1.96*sdcond,
         lo = mcond - 1.96*sdcond) %>% 
  ggplot() +
  geom_point(aes(x = condition, y = mcond), size = 3, color = "blue") +
  geom_hline(yintercept = 0) +
  geom_errorbar(aes(x = condition, ymin = lo, ymax = hi), color = "blue") +
  labs(x = "Condition", y = "ATE", title = "Weight of Advice") +
  coord_flip() + theme_bw()
m2sim95

# Combine graphs for Figure 1
grid.arrange(m1sim95, m2sim95, ncol = 1)

#### Create a table of the results from Figure 1
# Table A8 in SI
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r+1,nrow(existingDF)+1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}

# Create a baseline table
Tables <- stargazer(model1b, model2b, style="ajps", 
                    title="Models of Trust in Automation", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels=c( "Average of Humans", "Algorithm Aggregation", "Computer Algorithm")
)

# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)

# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 34 (after \hline after "constant")
r <- 17

# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & & \\\\"
hline <- "\\hline"
newline <- "\\\\"

# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r+1)
Tables <- insertrow(Tables,hline,r+2)
# Get number of unique values in each grouping
num.scenarios1 <- sapply(ranef(model1b),nrow)[1]
num.scenarios2 <- sapply(ranef(model2b),nrow)[1]

# Get standard deviation of the random effect
stddev.model1b.scenarios <- attributes(VarCorr(model1b)$scenario)$stddev
stddev.model2b.scenarios <- attributes(VarCorr(model2b)$scenario)$stddev

# Create a LaTex character row for the random effect
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, "&", num.scenarios2, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1b.scenarios, 3), "&", round(stddev.model2b.scenarios, 3), "\\\\")

# Add these lines to the table
Tables <- insertrow(Tables,number.of.scenarios,r+3)
Tables <- insertrow(Tables,stddev.scenarios,r+4)
Tables <- insertrow(Tables, hline, r+5)

# Write the table to a file that can be inserted into the document
write.table(Tables,file=(here("REtablepoliticalNC.tex")),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)


#### Models with alternative baseline condition
# Table A12 in SI
model1c <- lmer(distance ~ RandHumanTreat + AvgHumanTreat + LernerTreat + (1|scenario_num),
                data = experiments); summary(model1c)

model2c <- lmer(adviceWt ~ RandHumanTreat + AvgHumanTreat + LernerTreat + (1|scenario_num),
                data = experiments[experiments$adviceWt <= 1,]); summary(model2c)

model3c <- lmer(brier ~ RandHumanTreat + AvgHumanTreat + LernerTreat + (1|scenario_num),
                data = experiments); summary(model3c)

# Create a table of the results
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r+1,nrow(existingDF)+1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}

# Create a baseline table
Tables <- stargazer(model1c, model2c, model3c, style="ajps", 
                    title="Models of Trust in Automation", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels=c("Random Human", "Average of Humans", "Computer Algorithm")
)

# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)

# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 34 (after \hline after "constant")
r <- 17

# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & & \\\\"
hline <- "\\hline"
newline <- "\\\\"

# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r+1)
Tables <- insertrow(Tables,hline,r+2)
# Get number of unique values in each grouping
num.scenarios1 <- sapply(ranef(model1c),nrow)[1]
num.scenarios2 <- sapply(ranef(model2c),nrow)[1]
num.scenarios3 <- sapply(ranef(model3c),nrow)[1]

# Get standard deviation of the random effect
stddev.model1c.scenarios <- attributes(VarCorr(model1c)$scenario)$stddev
stddev.model2c.scenarios <- attributes(VarCorr(model2c)$scenario)$stddev
stddev.model3c.scenarios <- attributes(VarCorr(model3c)$scenario)$stddev

# Create a LaTex character row for the random effect
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, "&", num.scenarios2, "&", num.scenarios3, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1c.scenarios, 3), "&", round(stddev.model2c.scenarios, 3), "&", round(stddev.model3c.scenarios, 3), "\\\\")

# Add these lines to the table
Tables <- insertrow(Tables,number.of.scenarios,r+3)
Tables <- insertrow(Tables,stddev.scenarios,r+4)
Tables <- insertrow(Tables, hline, r+5)

# Write the table to a file that can be inserted into the document
write.table(Tables,file=(here("REtablepolAltB.tex")),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)


#### Models with additional control variables
# Table A.7.1 in the SI
model1 <- lmer(distance ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + Anchoring + 
                 tia + age + ed + female +
                 (1|scenario_num),
               data = experiments); summary(model1)

model2 <- lmer(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + 
                 tia + age + ed + female +
                 (1|scenario_num),
               data = experiments[experiments$adviceWt <= 1,]); summary(model2)

# Create Table A7.1.1
# Function for adding rows to a data frame (will be used later)
insertrow <- function(existingDF, newrow, r) {
  existingDF[seq(r+1,nrow(existingDF)+1),] <- existingDF[seq(r,nrow(existingDF)),]
  existingDF[r,] <- newrow
  existingDF
}

# Create a baseline table
Tables <- stargazer(model1, model2, style = "ajps", 
                    title = "Models of Trust in Automation", 
                    dep.var.labels.include = FALSE, 
                    covariate.labels=c( "Average Human Treatment", "Algorithm + Human Treatment", "Algorithm Treatment", "Anchoring Treatment","Trust in Automation","Age", "Education", "Female")
)

# Convert Tables to data frame and coerce the one vector (also called Tables) to character
Tables <- as.data.frame(Tables)
Tables$Tables <- as.character(Tables$Tables)

# Find where you want to put in the random effect. In our case, this is right after the last fixed effect. Line: 34 (after \hline after "constant")
r <- 27

# Create some standard label lines
randomeffect <- "{\\bf Random Effect} & & \\\\"
hline <- "\\hline"
newline <- "\\\\"

# Now insert those label lines where they are needed
Tables <- insertrow(Tables, hline, r)
Tables <- insertrow(Tables,randomeffect,r+1)
Tables <- insertrow(Tables,hline,r+2)
# Get number of unique values in each grouping
num.scenarios1 <- sapply(ranef(model1),nrow)[1]
num.scenarios2 <- sapply(ranef(model2),nrow)[1]
# Get standard deviation of the random effect
#stddev.model1.participants <- attributes(VarCorr(model1)$ResponseId)$stddev
stddev.model1.scenarios <- attributes(VarCorr(model1)$scenario)$stddev
#stddev.model2.participants <- attributes(VarCorr(model2)$ResponseId)$stddev
stddev.model2.scenarios <- attributes(VarCorr(model2)$scenario)$stddev

# Create a LaTex character row for the random effect
#number.of.participants <- paste("\\# of Participants & ", num.participants1, "&", num.participants2, "&", num.participants, "&", num.participants, "&", num.participants, "&", num.participants, "&", num.participants, "&", num.participants, "\\\\")
#stddev.participants <- paste("Participant Standard Deviation & ", round(stddev.model1.participants, 3), "&", round(stddev.model2.participants, 3), "&", round(stddev.model3.participants, 3), "&", round(stddev.model4.participants, 3), "&", round(stddev.model5.participants, 3), "&", round(stddev.model6.participants, 3), "&", round(stddev.model7.participants, 3), "&", round(stddev.model8.participants, 3), "\\\\")
number.of.scenarios <- paste("\\# of Scenarios & ", num.scenarios1, "&", num.scenarios2, "\\\\")
stddev.scenarios <- paste("Scenario Standard Deviation & ", round(stddev.model1.scenarios, 3), "&", round(stddev.model2.scenarios, 3), "\\\\")

# Add these lines to the table
#Tables <- insertrow(Tables,number.of.participants,r+3)
#Tables <- insertrow(Tables,stddev.participants,r+4)
Tables <- insertrow(Tables,newline,r+3)
Tables <- insertrow(Tables,number.of.scenarios,r+4)
Tables <- insertrow(Tables,stddev.scenarios,r+5)
Tables <- insertrow(Tables, hline, r+6)

# Write the table to a file that can be inserted into the document
write.table(Tables, file = ("REpol(wcov)gp.tex"),
            sep="",row.names= FALSE,na="", quote = FALSE, col.names = FALSE)


#### Heterogeneous treatment effects
# Left-hand side of Figure A22 in SI
# algorithm * tia INXN
model5 <- lmer(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + 
                 I(LernerTreat*tia) + 
                 tia + age + ed + female + 
                 (1|scenario_num), 
               data = experiments[experiments$adviceWt <= 1,]); summary(model5)

# algorithm * age INXN
model8 <- lmer(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + 
                 I(LernerTreat*age) +
                 tia + age + ed + female + 
                 (1|scenario_num), 
               data = experiments[experiments$adviceWt <= 1,]); summary(model8)

# algorithm * ed INXN
model11 <- lmer(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + 
                  I(LernerTreat*ed) + 
                  tia + age + ed + female + 
                  (1|scenario_num), 
                data = experiments[experiments$adviceWt <= 1,]); summary(model11)

# algorithm * female INXN
model14 <- lmer(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + 
                  I(LernerTreat*female) + 
                  tia + age + ed + female + 
                  (1|scenario_num), 
                data = experiments[experiments$adviceWt <= 1,]); summary(model14)

# INXN PLOTS:
# Model 5 Interaction Plot: TIA x Algorithm (weight)
mod5sum <- summary(model5)
tiaX <- seq(1,7)
maineffect5 <- fixef(model5)[4] + fixef(model5)[5] * tiaX
maineffect5 <- as.tibble(cbind(tiaX,maineffect5))
mod5sim <- sim(model5, n.sims = 1000)
mod5simcoef <- as.tibble(data.frame(coef(mod5sim)$fixef))
mod5simcoef$simulation <- rownames(mod5simcoef)
mod5simcoeflong <- mod5simcoef %>%
  mutate(simeffect1 = LernerTreat + I.LernerTreat...tia. * 1,
         simeffect2 = LernerTreat + I.LernerTreat...tia. * 2,
         simeffect3 = LernerTreat + I.LernerTreat...tia. * 3,
         simeffect4 = LernerTreat + I.LernerTreat...tia. * 4,
         simeffect5 = LernerTreat + I.LernerTreat...tia. * 5,
         simeffect6 = LernerTreat + I.LernerTreat...tia. * 6,
         simeffect7 = LernerTreat + I.LernerTreat...tia. * 7) %>%
  dplyr::select(simulation:simeffect7) %>%
  gather(tialevel, estimate, simeffect1:simeffect7) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel), nchar(tialevel))))

p4 = ggplot() + geom_line(data = mod5simcoeflong,aes(x=tia,y=estimate, group=simulation), color = "red", alpha = 0.1) +
  geom_line(data = maineffect5, aes(x=tiaX, y=maineffect5), size = 1.5) +
  xlab("Trust in Automation") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p4

# Model 8 Interaction Plot: Age x Algorithm (weight)
mod8sum <- summary(model8)
tiaX <- seq(19,74, by = 5)
maineffect8 <- fixef(model8)[4] + fixef(model8)[5] * tiaX
maineffect8 <- as.tibble(cbind(tiaX,maineffect8))
mod8sim <- sim(model8, n.sims = 1000)
mod8simcoef <- as.tibble(data.frame(coef(mod8sim)$fixef))
mod8simcoef$simulation <- rownames(mod8simcoef)
mod8simcoeflong <- mod8simcoef %>%
  mutate(simeffect19 = LernerTreat + I.LernerTreat...age. * 19,
         simeffect24 = LernerTreat + I.LernerTreat...age. * 24,
         simeffect29 = LernerTreat + I.LernerTreat...age. * 29,
         simeffect34 = LernerTreat + I.LernerTreat...age. * 34,
         simeffect39 = LernerTreat + I.LernerTreat...age. * 39,
         simeffect44 = LernerTreat + I.LernerTreat...age. * 44,
         simeffect49 = LernerTreat + I.LernerTreat...age. * 49,
         simeffect54 = LernerTreat + I.LernerTreat...age. * 54,
         simeffect59 = LernerTreat + I.LernerTreat...age. * 59,
         simeffect64 = LernerTreat + I.LernerTreat...age. * 64,
         simeffect69 = LernerTreat + I.LernerTreat...age. * 69,
         simeffect74 = LernerTreat + I.LernerTreat...age. * 74) %>%
  dplyr::select(simulation:simeffect74) %>%
  gather(tialevel, estimate, simeffect19:simeffect74) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel)-1, nchar(tialevel))))

p6 = ggplot() + geom_line(data = mod8simcoeflong,aes(x=tia,y=estimate, group=simulation), color = "red", alpha = 0.1) +
  geom_line(data = maineffect8, aes(x=tiaX, y=maineffect8), size = 1.5) +
  xlab("Age") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p6

# Interaction graph for model 11: Education x Algorithm (weight)
mod11sum <- summary(model11)
tiaX <- seq(2,4)
maineffect11 <- fixef(model11)[4] + fixef(model11)[5] * tiaX
maineffect11 <- as.tibble(cbind(tiaX,maineffect11))
mod11sim <- sim(model11, n.sims = 1000)
mod11simcoef <- as.tibble(data.frame(coef(mod11sim)$fixef))
mod11simcoef$simulation <- rownames(mod11simcoef)
mod11simcoeflong <- mod11simcoef %>%
  mutate(simeffect2 = LernerTreat + I.LernerTreat...ed. * 2,
         simeffect3 = LernerTreat + I.LernerTreat...ed. * 3,
         simeffect4 = LernerTreat + I.LernerTreat...ed. * 4) %>%
  dplyr::select(simulation:simeffect4) %>%
  gather(tialevel, estimate, simeffect2:simeffect4) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel), nchar(tialevel))))

p8 = ggplot() + geom_line(data = mod11simcoeflong,aes(x=tia,y=estimate, group=simulation), color = "red", alpha = 0.1) +
  geom_line(data = maineffect11, aes(x=tiaX, y=maineffect11), size = 1.5) +
  xlab("Education") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p8

# Interaction graph for model 14: Female x Algorithm (weight)
mod14sum <- summary(model14)
tiaX <- seq(1,2)
maineffect14 <- fixef(model14)[4] + fixef(model14)[5] * tiaX
maineffect14 <- as.tibble(cbind(tiaX,maineffect14))
mod14sim <- sim(model14, n.sims = 1000)
mod14simcoef <- as.tibble(data.frame(coef(mod14sim)$fixef))
mod14simcoef$simulation <- rownames(mod14simcoef)
mod14simcoeflong <- mod14simcoef %>%
  mutate(simeffect1 = LernerTreat + I.LernerTreat...female. * 1,
         simeffect2 = LernerTreat + I.LernerTreat...female. * 2) %>%
  dplyr::select(simulation:simeffect2) %>%
  gather(tialevel, estimate, simeffect1:simeffect2) %>%
  mutate(tia = as.numeric(substring(tialevel, nchar(tialevel), nchar(tialevel))))

p10 = ggplot() + geom_line(data = mod14simcoeflong,aes(x=tia,y=estimate, group=simulation), color = "red", alpha = 0.1) +
  geom_line(data = maineffect14, aes(x=tiaX, y=maineffect14), size = 1.5) +
  xlab("Female") + ylab("Estimated Impact of Algorithm") + ggtitle("Weight of Advice") +
  geom_hline(aes(yintercept = 0)) + theme_bw()
p10


#### Separate models for each of the geopolitical experiments
# Table A22 in SI
m1.w1.distance <- lm(distance ~ AvgHumanTreat + AlgHumanTreat + LernerTreat, 
                     data = experiments[experiments$scenario_num == 1,]); summary(m1.w1.distance)

m1.w2.distance <- lm(distance ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + Anchoring, 
                     data = experiments[experiments$scenario_num == 2,]); summary(m1.w2.distance)

m1.w3.distance <- lm(distance ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + Anchoring, 
                     data = experiments[experiments$scenario_num == 3,]); summary(m1.w3.distance)

m1.w4.distance <- lm(distance ~ AvgHumanTreat + AlgHumanTreat + LernerTreat + Anchoring, 
                     data = experiments[experiments$scenario_num == 4,]); summary(m1.w4.distance)

# Table A23 in SI
wave1 <- subset(experiments, scenario_num == 1)
wave2 <- subset(experiments, scenario_num == 2)
wave3 <- subset(experiments, scenario_num == 3)
wave4 <- subset(experiments, scenario_num == 4)

m1.w1.weight <- lm(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat, 
                   data = wave1[wave1$adviceWt <= 1,]); summary(m1.w1.weight)

m1.w2.weight <- lm(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat, 
                   data = wave2[wave2$adviceWt <= 1,]); summary(m1.w2.weight)

m1.w3.weight <- lm(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat, 
                   data = wave3[wave3$adviceWt <= 1,]); summary(m1.w3.weight)

m1.w4.weight <- lm(adviceWt ~ AvgHumanTreat + AlgHumanTreat + LernerTreat, 
                   data = wave4[wave4$adviceWt <= 1,]); summary(m1.w4.weight)


